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Abstract. At low reduced electric fields the electron energy distribution function in 
heavy noble gases can take two distinct shapes. This “bistability effect” - in which 
electron-electron (Coulomb) collisions play an essential role - is analyzed here for Xe 
with a Boltzmann equation approach and with a first principles particle simulation 
method. The solution of the Boltzmann equation adopts the usual approximations 
of (i) searching for the distribution function in the form of two terms (“two-term 
approximation”), (ii) neglecting the Coulomb part of the collision integral for the 
anisotropic part of the distribution function, (iii) treating Coulomb collisions as 
binary events, and (iv) truncating the range of the electron-electron interaction 
beyond a characteristic distance. The particle-based simulation method avoids these 
approximations: the many-body interactions within the electron gas with a true (un¬ 
truncated) Coulomb potential are described by a Molecular Dynamics algorithm, while 
the collisions between electrons and the background gas atoms are treated with Monte 
Carlo simulation. We find a good general agreement between the results of the two 
techniques, which confirms, to a certain extent, the approximations used in the solution 
of the Boltzmann equation. The differences observed between the results are believed 
to originate from these approximations and from the presence of statistical noise in 
the particle simulations. 
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1. Introduction 

The physics of electron swarms has been attracting considerable attention dnring the 
past decades becanse of the interest in the varions physical effects taking place in these 
settings and by the need for accnrate inpnt data in discharge modeling [HE]. The 
advance of the theoretical backgronnd, as well as of the compntational resonrces and 
nnmerical techniqnes made it possible to develop a detailed pictnre of the physics of 
particle swarms. Most of the efforts have been devoted to electron swarms, the electron 
energy distribntion fnnction and the transport properties have been determined for 
a wide variety of gases and gas mixtnres, and for a broad range of conditions, see 
e.g. [3HI3]. 

One of the particnlar phenomena, the bistability of the electron energy distribution 
function (EEDF) is a prononnced manifestation of the nonlinear natnre of these systems. 
The term “EEDF bistability” is nsed here to designate a sitnation when in a physical 
system, nnder hxed conditions (gas nnmber density, gas temperatnre, electric held 
strength, electron and ion concentrations, popnlation of excited states), two stable 
steady states are possible with different EEDFs. 

The electron temperatnre bistability was (according to onr best knowledge) predicted 
for the hrst time for a plasma nnder the effect of an alternating electromagnetic held, at 
conditions when electron-ion collisions are dominant mils]. Assnming a Maxwellian 
EEDF, the bistability in the plasma of heavy rare gases (argon, krypton and xenon) 
with an applied electric held, was shown to exist in ng dne to the specihc shapes 
of the momentnm transfer cross sections for the electron scattering from Ar, Kr and 
Xe atoms. Later, the ehect in heavy rare gases plasma was conhrmed by nsing a 
more accnrate approach, the Boltzmann eqnation (BE) analysis that took into acconnt 
electron-electron (e-e) collisions [IM9]. In nMn] it was fonnd that, within a certain 
range of parameters: (i) the rednced electric held strength, E/n and (ii) the electron to 
nentral atom density, n^jn (where E is the electric held strength, n is the atom nnmber 
density, and is the electron density), the BE has two stable solntions. In the case of 
xenon gas, e.g., the bistability ehect was fonnd at n^/n > 10“® and 0.025 Td < E/n < 
0.043 Td. 

The conditions mentioned above can be realized in diherent physical systems, e.g., in 
decaying plasmas in the presence of low electric held, in non-self-sustained discharges, 
and under electron swarm conditions. A decaying plasma in Xe was studied in Ref. in. 
where the time-dependence of the current in the afterglow plasma in the presence of a 
weak electric held was measured. A jump-like decrease in the current was observed at 
some time point during the plasma decay, which was attributed to the manifestation 
of the bistability ehect [16]. Non-self-sustained discharges (sustained by a beam of fast 
electrons) were considered in Refs. [20] 1^. in which the possibility of the existence of 
the EEDF bistability was theoretically analyzed in Ar, Kr and Xe. It was shown that 
the bistability ehect takes place in Xe and Kr, while for an Ar plasma [2T], the BE has 
a unique solution over the entire parameter range examined. To our best knowledge. 
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the bistability effect has not been investigated so far for electron swarm conditions. 

We note that besides pure noble gases, the bistability effect has also been studied in gas 
mixtures. In [22] the electron temperature (Te) in Ar:N 2 = 100:1 mixture afterglow was 
studied both experimentally and theoretically. In the experiments it was observed that, 
under certain conditions, a rather sharp knee appears in Te(t), while the vibrational 
temperature of nitrogen molecules, T^, remains almost constant. The theoretical study 
of the EEDF in Ar:N 2 afterglow plasma was carried out by the numerical solution of the 
appropriate BE, by taking into account e-e collisions, as well as superelastic vibrational 
and superelastic electronic collisions. The calculations have shown that ranges of rig and 
Tv exist, where two different solutions of BE can be obtained. The observed knee-like 
time-dependence of Tg was explained as the manifestation of the bistability effect. The 
possibility of this effect in pure nitrogen afterglow plasma was theoretically addressed 
in [23]. In addition, using the Maxwellian distribution function approach the bistability 
effect was studied in [23] for positron swarms in He. For more details the reader is 
referred to the review ESI- 

It should be noted that so far the theoretical analysis of the bistability effect was 
performed based on BE approach only. In these studies, and generally, in the 
solutions of the Boltzmann equation including electron-electron collisions, rather serious 
approximations have traditionally been used: the methods (i) search for the distribution 
function in the form of two terms (“two-term approximation”), (ii) neglect the electron- 
electron part of the collision integral for the anisotropic part of the distribution function, 
(iii) treat Coulomb collisions as binary events, and (iv) truncate the range of the electron- 
electron interaction beyond a characteristic distance. These approximations became 
widely accepted throughout the years, and one needs to note that in the absence of a 
more rigorous approach, their effects on the results of BE solutions have never been 
critically examined. To do that - and this is actually the motivation of our work - one 
would ideally describe the physical system of interest with a first-principles method that 
does not involve any approximations. (We note that previous particle-based (Monte 
Carlo) methods, dealing with Coulomb collisions, are neither free of approximations.) 
The particle-based simulation technique, used here [26], on the other hand, relaxes 
both the above major approximations and provides a solution based on Erst principles. 
Therefore we expect this method to be able to evaluate the accuracy of the BE solution 
that is limited by the approximations adopted, and can verify the effects predicted 
within the frame of the BE approach, like the EEDF bistability that was introduced 
above and is studied here in details in Xe gas. 

The physical system investigated here can be described by three independent 
parameters: the reduced electric held, E/n, the electron to gas (number) density ratio, 
nf,ln., and the electron density, rzg. In the description of the solution of the Boltzmann 
equation (section 3.2) this will be shown to follow from the structures of the equations 
involved. Equivalently, other sets of parameters derived from the above set can be used 
as well: we shall carry out our studies by varying E/n and Ue/n., while keeping the gas 
number density, n, constant. 
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In section 2 we introduce the underlying physical processes responsible for the 
development of the bistability effect. This is followed in section 3 by the presentation of 
the methods of calculation: section 3.1 discusses the particle simulation method, while 
section 3.2 outlines the solution of the Boltzmann equation. The results obtained by 
the two different approaches are presented and compared to each other in section 4. A 
short summary is given in section 5. 

2. Qualitative explanation of the EEDF bistability effect in rare gases 

The qualitative explanation of the possibility of the EEDF bistability effect in rare gases 
at weak electric helds (where only elastic collisions occur) adopts an approximation 
of a Maxwellian EEDF (i.e., it assumes a priori that the EEDF is Maxwellian, 
foie) ~ exp(—£//cBT'e), for details see e.g. [25]). For such conditions, the time-dependent 
equation for the electron temperature is as follows: 

irTi 

-^ = n <I>(E/n,Te) = n [H{E/n,T,) - L(Te,Tg)], (1) 

where the term H (Tg) describes the heating of electrons by the electric field and the term 
L(Te) accounts for the loss of electron energy in elastic collisions. (Note that, generally, 
the loss rate is positive, as electrons deposit energy in collisions with the background 
gas. However, at Tg < Tg, where Tg is the gas temperature, electrons are heated by the 
atoms of the gas, resulting in a negative loss, T < 0.) 

The number of steady state and stable solutions of ([ID depends on the shape of the 
<h(Tg) function. Figure [T](a) shows this function together with the H{Te) and L{Tff) 
functions [25], for Xe at E/n = 0.035 Td, and a gas temperature of Tg = 300 K. Let 
us note that L{Tf) is proportional to the elastic momentum transfer cross section crm(£) 
of electron scattering from Xe atoms, while H{Tf) varies proportionally to l/crm(£) (see 
section 3a). Due to the specihc shape of (Jm(£) (a deep Ramsauer-Townsend minimum 
at an energy of ~ 0.6 eV, see £gure[T](b)) the H{Tf) and L{Tf) curves exhibit a “wavy” 
behavior and cross each other at three points. As a consequence, the ^(Tg) function has 
an “inverse-N-type” shape and is equal to zero at three Tg values. Among these three 
steady state solutions of equation ([ID two (where 4’(Tg) has a negative slope) are stable, 
as indicated in figure [I](a). 

Although the EEDF in reality is generally not Maxwellian, but the e-e collisions drive the 
EEDF towards Maxwellian, therefore the above explanation holds at least qualitatively. 
We note that e-e collisions are essential in establishing the bistability effect, as the 
terms corresponding to these collisions introduce non-linearity into the BE. Without e-e 
collisions the BE is linear equation (in /o). As a linear equation, it has a unique solution 
at fixed parameters (taking into the account the proper normalization and the condition 
that /o —)■ 0 at electron energy oo). If e-e collisions are involved, the BE becomes a 
non-linear equation and, in principle, two or more different solutions are possible. For 
multiple solutions to appear the e-e collision frequency should be high enough to have 
an influence on the shape of the EEDF, i.e., it should be comparable or higher than the 
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frequency of energy losses in electron-atom elastic collisions. We note, however, that 
the non-linearity is a necessary, but not a sufficient condition for BE to have multiple 
solutions. While some examples of physical systems, where bistability is present, were 
given in section 1, a general sufficient condition for the appearance of the effect cannot 
be formulated at present. 



e[eV] 


Figure 1. (a) Heating (iJ) and loss (L) rates of the electron energy, due to the effect of 
the electric field and to their interaction with the background gas atoms, respectively. 
The curves have been derived assuming a Maxwellian EEDF and a gas temperature of 
Tg = 300 K. The data are given in arbitrary units, the grey dots indicate the stable 
equilibrium solutions, with $ = H — L = 0. (Note that a similar plot appeared in 
Ref. [25], as figure 27. In that figure the labelling of the H and L curves was exchanged 
as a misprint.) (b) Elastic momentum transfer cross section for electron - Xe atom 
collisions m- 


3. Methods 

Here we describe the basics of the two computational approaches, the particle simulation 
scheme is presented in section 3.1, while the solution of the Boltzmann equation is 
discussed in section 3.2. 

3.1. Particle simulation 

A swarm of electrons is simulated under the effect of a homogeneous electric held, in 
xenon gas, at different values of the number density ratio, rj = Uf-jn. At the low reduced 
helds considered only elastic collisions take place. As this ensures conservation of the 
number of particles, we consider a spatially homogeneous system. 

The simulation scheme is based on a combination of a Molecular Dynamics (MD) 
technique and a Monte Carlo (MC) approach [26]. The MD describes the many-body 
interactions (driven by the inter-particle Coulomb potential) within the classical electron 
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gas, while the MC part handles the interaction of the electron gas with the backgronnd 
(atomic) gas. 


3.1.1. Molecular Dynamics simulation of the electron gas. In our classical many-body 
system the particles move under the influence of interparticle and external forces in a 
cubic simulation box, with periodic boundary conditions. The equations of motion of 
particles i = 1,... ,N, which are numerically integrated with discrete time steps (At), 
are Newtonian: 


m 


dD 


Fjj eE, 

Do 


( 2 ) 


where m and —e are the mass and the charge of the electron, respectively. The sum 
on the right hand side represents the force exerted on particle i by all other {j ^ i) 
particles and their periodic images located in spatial replicas of the simulation box. 
These images have to be included in the proper determination of the interparticle forces 
in the case of the “full” (un-truncated) inhnite-range Coulomb potential. (Note that for 
our conditions no screening of the potential takes place by oppositely charged species.) 
This summation is indeed a key issue and needs a special approach, as will be explained 
below. The second term on the right hand side of the above equation is a contribution 
due to the external electric held. We note that in the absence of the interaction of the 
electrons with a background gas the electrons would continuously be accelerated by this 
held, which will, however, not be the case when e“+Xe atom collisions take place. 

The MD method adopted in this study is based on the Particle-Particle Particle-Mesh 
(PPPM) approach, described in details in |2H], that uses a partitioning of the interaction 
into (i) a force component that can be calculated on a mesh (the “mesh force”) and (ii) 
a short-range (“correction”) force, which is to be applied to closely separated pairs of 
particles only. In the mesh part of the calculation, charge clouds are used instead of 
point-like charges. The charge density distribution p(r) is assigned to a grid and is 
Fourier transformed to the k-space. Multiplying p(k) with an optimized Green function 
results in a potential distribution 0(k) = G(k)p(k), which is subsequently transformed 
back to real space. The forces acting on the particles are obtained by differentiation 
of the potential and interpolating the electric held to the positions of the particles. 
The cloud shape is chosen in a way to ensure that p(k) is band-limited (which would 
not hold for point-like charges). The calculation of the potential in the Fourier space 
automatically takes into account the periodic images of the primary computational 
box. For closely separated particles a correction force is to be applied, which is the 
difference between the forces between two point-like charges and two charges with the 
cloud shape used. The simulations describe a micro-canonical ensemble, where the 
number of particles, the volume of the system and energy are conserved. 

The upper limit for the simulation time step is defined by the stability of the integration 
of the equation of motion in the case of the closest approach of two electrons, rmin = 
e^/ (dTreo^max)- Here ^max is a pre-dehned maximum energy [28], which has to be chosen 
carefully, to ensure that the probability of finding electrons with e > is vanishingly 
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small at the conditions considered. In this work we adopt a maximum energy of 5 eV, 
for this value the time step has to be chosen to be as low as At = 2.2 x 10“^® s. This 
results in a very demanding computational load to follow the evolution of the system 
for a long time. As a converged solution assumes sufficient interaction between the 
electrons and the gas as well, one needs to have a high-enough electron-atom collision 
frequency. This can only be ensured by setting a high the gas atom number density. 
For this reason our computations are carried out with n = 2.5 x 10^^ m“^, which allows 
simulated times of tens of nanoseconds (with a computational speed of simulating ~ 1 
ns in one day, on a single CPU). The number of electrons is chosen to be = 1000. 

3.1.2. Monte Carlo simulation of electron gas - background gas interaction. Having 
solved the description of the electron gas with the MD method that accounts for electron- 
electron interactions, now we introduce a background gas and let the two gases to 
interact via e“-|-Xe collisions. The probability of an e“-|-Xe collision during a time step 
At is calculated as: 

^coii = 1 - exp[-?r(Tm(fi')5'At], (3) 

where g = |g|, with g = v —V being the relative velocity between the electron and a Xe 
atom with a velocity V randomly chosen from a Maxwellian background of gas atoms 
having a temperature Tg. This probability is calculated for each electron in each time 
step, and decision about the occurrence of a collision is made by comparing Pcoii with a 
random number. Collisions are executed in the center-of-mass frame, and are considered 
to be isotropic. The energy change of the gas atoms colliding with the electrons is not 
accounted for, the gas temperature is kept constant. 

The simulations allow investigation of both the stationary state and the transients 
(induced by changing the electric held strength or the gas temperature). The transients 
are followed by monitoring the time-dependence of the mean energy of the electrons, 
which is calculated in each time step. The EEDF is, on the other hand, only calculated 
for the stationary state as its accurate determination (over 6-7 orders of magnitude) 
requires averaging over millions of time steps. 


3.2. Solution of the Boltzmann equation 


Here we discuss the specihc features of the BE for electrons, as applied to conditions 
under consideration (homogeneous plasma, or an electron swarm in an atomic gas, acted 
upon by a weak steady electric held). The distribution function of the electrons, /(v), 
can be described by the equation 
dfiv) eE 


Vv/(v) = U, 


( 4 ) 


dt m 

where v is the electron velocity and C is the collision integral. Further we shall 
consider only elastic scattering of electrons from atoms and electron-electron collisions: 

c = c^ + a. 
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The solution of the BE is based on expansion of the distribution function in Legendre 
polynomials P„(cos0), in which only two first terms are taken into account: 

/(v) =/o(n)+ /i(n)cos0, (5) 

where v is the velocity magnitude, 0 is the angle between v and —E, fo{v) is the 
symmetrical part of the distribution function and fi{v) describes the directed motion 
of the electrons along the electric field. (We note that this “two-term approximation”, 
in the absence of e-e collisions, has been benchmarked with other solution methods of 
the BE, as well as with particle-based (Monte Carlo) simulations in several studies, see 


The substitution of expansion ([5]) into equation (jlj) leads to equations for the /o and /i 
functions: 


and 


dfo 

dt 

dt 


eE d 
3mv‘^ dv 

eEdfo 


(v^fi) = Com + C 


Oe 


= + C: 


le* 


m dv 

The collision integrals Com and Cim can be written as |15j : 


Com — 


1 d 


2^2 dv 
Clm ~ ^mfli 


2m 

~M 


Z/mU 


ksT^dfo 
m dv 


+ vfo 


( 6 ) 

( 7 ) 

( 8 ) 
( 9 ) 


where z/m = na^v is the momentum transfer frequency and 6 = 2m/M is the average 
fraction of the energy lost by the electrons in one elastic collision with atom (M is the 
mass of the gas atom). The rate of the electron energy loss due to elastic collisions is 
characterized by the frequency = (5z/m- 

It is known that the calculation of the pair-collision frequency in the case of Coulomb 
collisions encounters a characteristic difficulty, namely the logarithmic divergence of 
frequency at small scattering angles. This difficulty is avoided by assuming that the 
Coulomb potential acts only up to a certain finite distance rmax (see later). The 
expression for the term S'oe is written as follows [T^[30] : 


^0e = 


dv 




Ai{fo)vfo + ^2(/o 


9fo 

dv 


Aiifo) =471 f vlfo{vi)dvi, 

JO 

rv roo 

/ v^fo{vi)dvi+v^ vifo{vi)dvi 
Jo Jv 


dvr 

Mfo) = y 


where fo{v) is normalized as: 

poo 

dvr / v‘^fo{v)dv = 1 
Jo 


and 


z/p = 27mp 




, ro = 


d'KEomv'^ 


( 10 ) 

( 11 ) 

( 12 ) 

(13) 

(14) 
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In BE calculations the value of the parameter tq is usually estimated by means of the 
mean electron energy, that gives cq = /‘iT^Sole . For the case of plasmas the Debye 
length, Ad, is taken to be the cutoff distance, i.e. rmax = Ad- For the case of swarm 
conditions considered here we use the approximation: 

?^max = (15) 


i.e. the half of the average distance between the electrons is taken to be the cutoff 
distance. (This choice of r^ax is based on intuitive physical considerations: if the impact 
parameter (of test electron relative to a given electron) is higher than the half the average 
distance between electrons in the gas, then the influence of this electron on the test one 
becomes weaker than the influence of the neighboring one.) Note that in most cases 
S> 1 (see comments in [TM30] i and the “1” in the expression under the logarithm 
in (ITT)) can be omitted. 

As to the term S'le in equation (4), it is very complex (see comments in [15]) and we 
did not hnd publications in which this term was taken into account in calculations. As 
a rule (BOLSIG+ [31], EEDF [32l[33]), it is neglected assuming that Ug Then, if 

the characteristic time of plasma parameters variation is essentially smaller than z/“^, 
the time derivative in equation (I?]) can be omitted. In this case 
eE dfo 


fi = 


z/mm dv 


and equation for the /o function is written as 


dfo 

dt 


eE d ( 2 eE dfo' 


dv\ Vmm dv 


— Gom + 


Oe- 


(16) 


( 17 ) 


It should be noted that, in the absence of e-e collisions, the parameter for the steady state 
solution of equation (ITT)) is the reduced electric held E/n. If the e-e collisions are taken 
into account, as it has already been mentioned in section 1, there are three parameters: 
E/n, nf,ln and Uf.. The electron number density is an independent parameter since the 
logarithmic term in eq. (ITT)) (the Coulomb logarithm) depends on rig. Actually, at hxed 
E/n and n^/n values the dependence of /o on n^ is rather weak, since n^ value is under 
logarithm. 

For the numerical solution of (1161) it is rewritten with energy as a variable. The steady- 
state equation is solved by an iteration method similar to that in [33l[3l]. The initial 
/o(e) is assumed to be the Maxwellian with a given electron temperature T^o- In the 
case of calculation of the time-dependent solution of the BE the time step should be as 
small as At -C and At For conditions under consideration At was taken as 

At = 2 X 10-1^ s. 


4. Results 

Bistability of the EEDF, i.e. two solutions for the distribution function, for certain 
conditions are found here, both via the solution of the Boltzmann equation and via 
executing the particle-based simulations. To hnd these different solutions both methods 
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Figure 2. (a) Mean energy values obtained from the calculations at different values 
of the electron to gas density ratio, g. The results of both the Boltzmann equation 
solutions (“BE”, lines) and of particle simulations (“PS”, symbols) originate from 
calculations started with different initial conditions. Tg = 300 K. (b) The effect of the 
gas temperature (Tg) on the limits of the bistable solutions, at 77 = 10“^. 


start from different initial conditions: both assnme a Maxwellian distribntion with a 
“low” and a “high” mean energy. For the particle method ^init = 0.009 eV and 0.9 
eV are nsed as “low-energy” and “high-energy” initial conditions, with random initial 
positions of the particles in the simnlation box. For the BE solntion the initial mean 
energy valnes are: Sinit = 0.065 eV and 0.65 eV. In the domain where a nniqne solntion 
exists the resnlts of the calculations do not depend on these choices. In the case of two 
possible solutions two domains of £init exist: from one of these domains (sinit < ^fnit) 
the calculations converge to the lower energy solution, and from the other domain 
(Sinit > ^hiit) to the higher energy solution. The value of Sinit within these domains 
has no effect on the results. The boundary value, depends on the E/n and 77 . 

The BE analysis has been carried out for a range of reduced electric helds, 0.015 Td 
< E/n < 0.05 Td, for various electron to gas density ratios. Results of these calculations, 
in terms of the mean electron energy e, for rj = 10“®, 10“^, and 10“® are displayed in 
£gure[2](a). Two solutions are found for the intermediate values of E/n, the boundaries 
of the domain change slightly with the density ratio. The “low-energy” solution exhibits 
a mean energy that is almost independent of r], and is in the range £ ~ 0.04... 0.05 eV. 
In the case of the “high-energy” solutions the mean energy amounts several tenth of 
an electron Volt and increases with increasing rj. The particle simulation results agree 
generally well with those obtained from the BE, small differences are found in the case 
of “high-energy” solutions. Additionally, the particle simulation does not predict a 
bistability at 77 = 10“® at E/n= 0.035 Td, whereas the BE solution does. The reason 
of this discrepancy is not fully understood. We contemplate that the presence of noise 
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Figure 3. (a,b,c) The convergence of the particle simulation method starting from 
a “high-energy” distribution (solid red lines) and a “low-energy” distribution (dashed 
dark blue lines), for the E/n values indicated. Two solutions are found for E/n = 
0.035 Td. The dotted horizontal black lines indicate the mean energy values obtained 
from BE calculations. (d,e,f) /(e) for the stable solutions: Boltzmann equation data 
(lines) in comparison with particle simulation data (symbols), y = 10”^ for all plots. 


in the particle simulation, due to its statistical nature, prevents finding solutions where 
the minimum in the energy balance vs. Tg is very shallow. This, however, needs further 
clarification. While all the data in figure [2](a) have been shown for a gas temperature of 
Tg = 300 K, figure |2](c) shows the dependence of the results on Tg, for fixed rj = 10“^. 
The mean energy here is displayed on a log scale, as the gas temperature influences 
predominantly the “low-energy” solutions. For a lower Tg a lower e is found, however, 
the domain of E/n, where bistability is found, is wider at Tg = 0 K. 

The stationary solutions in the particle simulations were typically obtained beyond a 
few tens of nanoseconds of simulated time. The convergence of the mean energy is 
illustrated in figures |3]( a)-(c), for reduced electric fields of E/n = 0.025 Td, 0.035 Td, 
and 0.045 Td, respectively, aX rj = 10“^. In all cases the simulations were started from 
two initial configurations, as already mentioned above. While at E/n = 0.025 Td and 
0.045 Td we observe convergence to a unique value of e, the runs at E/n = 0.035 Td 
clearly yield two solutions, with mean electron energies differing by almost an order of 
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magnitude. The full EEDFs, obtained by the two methods, are compared in hgures 
[3jd)-(f). The EEDFs are normalized as /q°° f{e)^de = 1. We hnd a good agreement 
between the data obtained from the BE and from the particle simulation, although 
some differences show up in the tails of the distribution functions. The difference of the 
EEDFs belonging to the two stable solutions at E/n = 0.035 Td is remarkable. 



f[ns] 

Figure 4. Time dependence of the mean electron energy (response of the swarm) 
following a sudden change of the electric field (at t = 0) from E/n = 0.035 Td to E'/n 
= 0.025 Td and 0.045 Td. The thick lines correspond to solutions of the BE, while thin 
lines show the particle simulation results. Unique solutions are found for the “new” 
values of the reduced electric field, y = 10“® 

To demonstrate the ability of our methods to follow the temporal evolution of the 
system here we take as an example the change of the electric held. We start from the 
two stationary stable solutions obtained at E/n = 0.035 Td (at 7] = 10“®), change 
E/n to a “new” value, and monitor the convergence of e{t). We test “new” values of 
E/n = 0.025 Td and 0.045 Td. The results obtained by both methods are depicted in 
m Following the change of the reduced held both methods - although with a diherent 
dynamics - show the relaxation of the system to unique states, corresponding to the 
new values of E/n. The relaxation time is found to be of the order of few times ten 
nanoseconds. Here we recall the importance of the energy loss frequency z/u, as this 
dehnes the time scale of the relaxation of the mean energy. An estimation at e = 0.5 
eV gives k, 20 ns, in accordance with our observations of the relaxation time scale. 
The peculiarities of the transitions can be understood by following the heating and 
the cooling of the electrons at the specihc conditions. The rate of electron heating 
is proportional to E"^/a^{e) and the rate of energy loss is proportional to (Ju(£) £ = 
{2m/M)ani{e) e. 

First we discuss the transition from the “low-energy” solution at E/n = 0.035 Td to the 
hnal state at E/n = 0.045 Td. The change of the mean energy in this case is slow at 
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the beginning but becomes very fast afterwards. In the initial state the EEDF is narrow 
and the mean energy is low, where the momentum transfer cross section is high. As the 
heating rate is proportional to and the cooling rate is proportional to crm(e), the 

rate of change of e is low at high in the first phase of the transition. When the mean 
energy becomes higher aya{e) decreases, giving more preference to the heating. This, 
together with the increasing slope of cr^is), results in an abrupt increase of e beyond a 
certain time. 

The transition from the “high-energy” solution at E/n = 0.035 Td to the final state at 
E/n = 0.025 Td is more complex. The three stages of relaxation - seen in figure H]- can 
be explained as follows. Under the steady state conditions at E/n = 0.035 Td the mean 
electron energy is relatively high, e k, 0.48 eV. At this energy the momentum transfer 
cross section is low (due to the Ramsauer-Townsend minimum) and, as a consequence, 
the rate of electron heating is relatively high. The high rate of heating is balanced by 
a high rate of energy loss. The decrease of the electric held causes an instant decrease 
of the heating. The unbalanced high rate of energy loss leads to the rapid decrease in 
the mean electron energy (during the hrst stage of e{t) relaxation) and, consequently, 
to a decrease of the rate of energy loss. The decrease in the rate of energy loss leads 
to near-balanced conditions, i.e. the difference between the rate of energy loss and the 
rate of heating becomes moderate, leading to the battening of the e{t) curve (second 
stage of relaxation). When the mean energy decreases below e ~ 0.2 eV, the momentum 
transfer cross section increases sharply and the rate of energy loss increases noticeably, 
while the rate of heating decreases signihcantly. Consequently, £ decreases faster below 
0.2 eV, during the third stage of the relaxation. We note that since under the conditions 
considered here the electron-electron collision frequency is higher than the electron-atom 
energy exchange frequency (r'u), the EEDF is nearly Maxwellian during the relaxation 
process. 

The transitions, which involve a small change of the mean energy (from the “high- 
energy” solution at E/n = 0.035 Td to the hnal state at E/n = 0.045 Td, and from 
the “low-energy” solution aX E/n = 0.035 Td to the hnal state at E/n = 0.025 Td) are 
much faster, compered to the two cases discussed above. 

5. Summary 

We have investigated the bistability of the EEDF in Xe gas at low reduced electric helds 
via the solution of the Boltzmann equation and via a hrst principles particle simulation 
technique. The solution of the Boltzmann equation adopted the usual, widely accepted 
approximations: it (i) searched for the distribution function in the form of two terms, 
(ii) neglected the electron-electron part of the collision integral for the anisotropic part 
of the distribution function, (iii) treated Coulomb collisions as binary events, and (iv) 
truncated the range of the electron-electron interaction beyond a characteristic distance. 
The particle simulation method [26], being devoid of any of these approximation has 
provided hrst-principles solutions to the problem, via a combination of a Molecular 
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Dynamics simulation method (that described accurately the many-body interactions 
within the electron gas governed by the full Coulomb potential) and a Monte Carlo 
method (that handled the interaction of the electrons with the atoms of the background 
gas). Both methods allowed the computation of the EEDF and the related quantities, 
and have indicated the existence of two stable solutions for the EEDF for a range of 
E/n. The electron mean energies and the full EEDFs, obtained by the two methods 
agreed generally well for most od the parameter settings covered. Differences found for 
the domain of the bistability and for the shapes of the EEDFs may be attributed to 
approximations adopted in the BE solutions and to the presence of statistical noise in 
the particle simulations. 
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